<!DOCTYPE html>
<html>
<head>
  <meta charset="utf-8">
  
  <title>OpenFOAM 中求解 ODE 一例 | Giskard&#39;s CFD Learning Tricks</title>
  <meta name="viewport" content="width=device-width, initial-scale=1, maximum-scale=1">
  <meta name="description" content="本篇介绍如何编写一个小程序来调用 OpenFOAM 的 ODE 求解器来求解任意常微分方程的初值问题。">
<meta property="og:type" content="article">
<meta property="og:title" content="OpenFOAM 中求解 ODE 一例">
<meta property="og:url" content="http://xiaopingqiu.github.io/2016/08/21/ODE/index.html">
<meta property="og:site_name" content="Giskard's CFD Learning Tricks">
<meta property="og:description" content="本篇介绍如何编写一个小程序来调用 OpenFOAM 的 ODE 求解器来求解任意常微分方程的初值问题。">
<meta property="og:image" content="/image/RKCK45.png">
<meta name="twitter:card" content="summary">
<meta name="twitter:title" content="OpenFOAM 中求解 ODE 一例">
<meta name="twitter:description" content="本篇介绍如何编写一个小程序来调用 OpenFOAM 的 ODE 求解器来求解任意常微分方程的初值问题。">
  
  
    <link rel="icon" href="/favicon.png">
  
  <link href="//fonts.googleapis.com/css?family=Source+Code+Pro" rel="stylesheet" type="text/css">
  <link rel="stylesheet" href="/css/style.css" type="text/css">
  
<!-- Google Analytics -->
<script type="text/javascript">
(function(i,s,o,g,r,a,m){i['GoogleAnalyticsObject']=r;i[r]=i[r]||function(){
(i[r].q=i[r].q||[]).push(arguments)},i[r].l=1*new Date();a=s.createElement(o),
m=s.getElementsByTagName(o)[0];a.async=1;a.src=g;m.parentNode.insertBefore(a,m)
})(window,document,'script','//www.google-analytics.com/analytics.js','ga');

##ga('create', '[object Object]', 'auto');
ga('create', 'UA-62501686-1', 'auto');
ga('send', 'pageview');

</script>
<!-- End Google Analytics -->


</head>
<body>
  <div id="container">
    <div id="wrap">
      <header id="header">
  <div id="banner"></div>
  <div id="header-outer" class="outer">
    <div id="header-title" class="inner">
      <h1 id="logo-wrap">
        <a href="/" id="logo">Giskard&#39;s CFD Learning Tricks</a>
      </h1>
      
        <h2 id="subtitle-wrap">
          <a href="/" id="subtitle">CFD and Scientific Computing</a>
        </h2>
      
    </div>
    <div id="header-inner" class="inner">
      <nav id="main-nav">
        <a id="main-nav-toggle" class="nav-icon"></a>
        
          <a class="main-nav-link" href="/">Home</a>
        
          <a class="main-nav-link" href="/archives">Archives</a>
        
          <a class="main-nav-link" href="/atom.xml">Rss</a>
        
          <a class="main-nav-link" href="/about">About</a>
        
      </nav>
      <nav id="sub-nav">
        
        <a id="nav-search-btn" class="nav-icon" title="Search"></a>
      </nav>
      <div id="search-form-wrap">
        <form action="//google.com/search" method="get" accept-charset="UTF-8" class="search-form"><input type="search" name="q" results="0" class="search-form-input" placeholder="Search"><button type="submit" class="search-form-submit">&#xF002;</button><input type="hidden" name="q" value="site:http://xiaopingqiu.github.io"></form>
      </div>
    </div>
  </div>
</header>

      <div class="outer">
        <section id="main"><article id="post-ODE" class="article article-type-post" itemscope itemprop="blogPost">
  <div class="article-meta">
    <a href="/2016/08/21/ODE/" class="article-date">
  <time datetime="2016-08-21T07:26:43.000Z" itemprop="datePublished">2016-08-21</time>
</a>
    
  <div class="article-category">
    <a class="article-category-link" href="/categories/OpenFOAM/">OpenFOAM</a>
  </div>

  </div>
  <div class="article-inner">
    
    
      <header class="article-header">
        
  
    <h1 class="article-title" itemprop="name">
      OpenFOAM 中求解 ODE 一例
    </h1>
  

      </header>
    
    <div class="article-entry" itemprop="articleBody">
      
        <p>本篇介绍如何编写一个小程序来调用 OpenFOAM 的 ODE 求解器来求解任意常微分方程的初值问题。 </p>
<a id="more"></a>
<h5 id="1-_数学背景">1. 数学背景</h5><p>首先简要看一下涉及的数学背景。对于一阶的常微分方程，<br>$$<br>y’=f(x,y), \quad x\in[a,b] \\<br>y(a)=y_0<br>$$<br>常微分方程，如果存在解析解的话，其解应该是一个函数 $y=f(x)$。然而，大多数常微分方程是没有解析解的，只能数值求解。数值方法得到的，是一系列的 $x_0, x_1, \cdots x_n$ 对应的函数值 $y_0, y_1, \cdots y_n$。</p>
<p>常用的数值解法有：</p>
<ul>
<li>显式欧拉法<br>这种方法最简单，将区间 $[a,b]$ 分成 n 份，则得到步长 $h=(b-a)/n$。以 $y(a)=y_0$ 为起点，通过下述迭代，<br>$$<br>y_{m+1} = y_m+hf(x_m,y_m)<br>$$<br>可以得到 $x_m=a+m*h$ 处的函数值 $y_m$，此即为常微分方程的数值解。显式欧拉法形式简单，但是只有一阶精度，而且稳定性是有条件的，一般在实际中较少用到。</li>
<li>改进的欧拉法<br>这种方法是在显式欧拉法的基础上改进得到，将显式欧拉法中使用的向前积分改为梯形积分。其迭代形式为<br>$$<br>y_{m+1} = y_m +\frac{h}{2}[f(x_m,y_m)+f(x_{m+1},y_{m+1})]<br>$$<br>这种方法不是显式地，所以，需要在每一步内进行迭代求解。可以用如下的迭代公式<br>$$<br>y_{m+1}^{(n+1)} = y_m +\frac{h}{2}[f(x_m,y_m)+f(x_{m+1},y_{m+1}^{(n)})]<br>$$<br>注意，这里的 n 指的是计算 $y_{m+1}$ 的值时的内迭代次数，迭代的初始值 $y_{m+1}^{(0)}$ 可以用显式欧拉公式来给出<br>$$<br>y_{m+1}^{(0)} = y_m + hf(x_m, y_m)<br>$$</li>
<li>显式 4 阶 Runge-Kutta 方法<br>这是现实中常用的一种方法，其迭代形式如下<br>$$<br>\begin{align*}<br>y_{m+1} &amp; = y_m + \frac{h}{6}[k_1+2k_2+2k_3+k_4] \\<br>k_1 &amp; = f(x_m, y_m) \\<br>k_2 &amp; = f(x_m+\frac{1}{2}h, y_m+\frac{1}{2}hk_1) \\<br>k_3 &amp; = f(x_m+\frac{1}{2}h, y_m+\frac{1}{2}hk_2) \\<br>k_4&amp;  = f(x_m+h, y_m+hk_3)<br>\end{align*}<br>$$</li>
</ul>
<p>除了以上，当然还有很多方法，比如预测校正等等，这里就不再逐一介绍了。</p>
<p>问题是，实际中遇到的还可能是高阶的常微分方程，比如，弹簧-谐振子系统可以用以下e二阶常微分方程描述<br>$$<br>\frac{d^2y}{dt^2} = -\frac{k}{m}y<br>$$</p>
<p>高阶常微分方程的初值问题，可以用以下通式来描述<br>$$<br>y^{n} = f(x,y,y’,\cdots y^{n-1})<br>$$</p>
<p>其初始条件为<br>$$<br>y(0) = a_0, \quad y’(0) = a_1, \quad \cdots, \quad y^{n-1}(0) = a_n<br>$$<br>对于这种高阶常微分方程，可以将其表述为一系列一阶常微分方程的组成的方程组来求解。下面以二阶常微分方程为例，介绍如何将高阶常微分方程的初值问题转化为一阶常微分方程组。<br>考察如下二阶常微分方程<br>$$<br>y’’ = f(x,y,y’), \quad x\in[a,b] \\<br>y(a) = a_0, \quad y’(a) = a_1<br>$$</p>
<p>若令 $z=y’$，则上述二阶常微分方程可以表示成如下方程组<br>$$<br>\left \{<br>\begin{align*}<br>y’ &amp;= z \\<br>z’ &amp;= f(x,y,z)<br>\end{align*}<br>\right. \\<br>y(a)=a_0,\quad z(a)=a_1<br>$$</p>
<p>这个方程组，就可以用前面介绍的一阶常微分方程的解法来求解了，比如，若用最简单的显式欧拉法，则<br>$$<br>\begin{align*}<br>y_{m+1} &amp;= y_m + h z_m \\<br>z_{m+1} &amp;= z_m +hf(x_m, y_m, z_m)<br>\end{align*}<br>$$<br>或者用显式 4 阶 Runge-Kutta 方法<br>$$<br>\begin{align*}<br>y_{m+1} &amp; = y_m + \frac{h}{6}[K_1+2K_2+2K_3+K_4] \\<br>z_{m+1} &amp; = z_m + \frac{h}{6}[M_1+2M_2+2M_3+K_4] \\<br>K_1 &amp; = z_m,\quad M_1=f(x_m, y_m, z_m) \\<br>K_2 &amp; = z_m + \frac{M_1}{2}, \quad M_2 = f(x_m+\frac{h}{2}, y_m+\frac{K_1}{2}, z_m+\frac{M_1}{2})\\<br>K_3 &amp; = z_m + \frac{M_2}{2},\quad M_3 = f(x_m+\frac{h}{2}, y_m+\frac{K_2}{2}, z_m+\frac{M_2}{2})\\<br>K_4&amp;  = z_m + M_3,\quad M_4=f(x_m+h, y_m+K_3, z_m+M_3)<br>\end{align*}<br>$$</p>
<p>二阶以上的常微分方程，除了可以给出初值条件，还可以给出边值条件，比如<br>$$<br>y’’ = f(x,y,y’), \quad x\in[a,b] \\<br>y(a) = \alpha, \quad y(b) = \beta<br>$$</p>
<p>这种情况下，就无法直接将此方程转化为一阶常微分方程组了。但是，边值问题可以通过一定的方法转换成初值问题，以下给出一种：<strong>试射法</strong>。<br>在不知道 $y’(a)$ 的情况下，不妨假设 $y’(a)=\gamma_1$，这样，就得到了一个初值问题<br>$$<br>y’’ = f(x,y,y’), \quad x\in[a,b] \\<br>y(a) = \alpha, \quad y’(a) = \gamma_1<br>$$<br>解此初值问题，得到 $y(b)$ 的值 $\beta_1$，并与 $\beta$ 比较，如果误差足够小，则认为假设的 $y’(a)=\gamma_1$ 是合理的。否则，就对 $\gamma_1$ 进行修正，比如令 $\gamma_2 = \tfrac{\beta}{\beta_1}\gamma_1$，然后再以 $y’(a)=\gamma_2$ 为初值，继续求解初值问题。直到得到的初值问题的解$y(b)=\beta_k$ 与 $\beta$ 足够接近为止。</p>
<p>上述方程可以归纳为，将初值问题转化为如下边值问题<br>$$<br>y’’ = f(x,y,y’), \quad x\in[a,b] \\<br>y(a) = \alpha, \quad y’(a) = \gamma_k, k=1,2,\cdots<br>$$<br>若记问题 $y_k(x)$ 的解为 $y(x;\gamma_k)$，则 $\gamma_k$ 的理想值应该满足<br>$$<br>F(\gamma) = y(b;\gamma)-\beta = 0<br>$$<br>这个方程，可以用牛顿迭代法来求解：<br>$$<br>\gamma_{k+1} = \gamma_k-\frac{F(\gamma_k)}{F’(\gamma_k)}<br>$$<br>其中，$F(\gamma_k)=y(b;\gamma_k)-\beta=\beta_k-\beta$。那么 $F’(\gamma_k)$ 该如何得到呢？根据 $F(\gamma_k)$ 的定义，可以知道 $F’(\gamma_k) = \frac{\partial y(b;\gamma)}{\partial \gamma}\big|_{\gamma=\gamma_k}$，若定义 $W=\frac{\partial y(b;\gamma)}{\partial \gamma}$，则 $F’(\gamma_k)=W(b;\gamma_k)$。</p>
<p>将上述归纳形式的初值问题，对$\gamma$ 求偏导，得<br>$$<br>\frac{\partial y’’}{\partial \gamma} = \frac{\partial f(x,y(x;\gamma),y’(x,y’))}{\partial y} \frac{\partial y(x;\gamma)}{\partial \gamma} + \frac{\partial f(x,y(x;\gamma),y’(x,y’))}{\partial y’} \frac{\partial y’(x,y’)}{\partial \gamma}<br>$$</p>
<p>根据 $W$ 的定义，有<br>$$<br>W=\frac{\partial y(x;\gamma)}{\partial \gamma}, W’=\frac{\partial y’(x,y’)}{\partial \gamma}, W’’=\frac{\partial y’’}{\partial \gamma}<br>$$</p>
<p>于是，可以得到一个关于 $W$ 的二阶常微分方程<br>$$<br>W’’=\frac{\partial f(x,y,y’)}{\partial y} W + \frac{\partial f(x,y,y’)}{\partial y’}W’<br>$$<br>其定解条件为<br>$$<br>W(a)=\frac{\partial y(a;\gamma)}{\partial \gamma}=0, W’(a)=\frac{\partial y’(a;\gamma)}{\partial \gamma}=\frac{\partial \gamma}{\partial \gamma} = 1<br>$$</p>
<p>这样就构成了一个关于 $W$ 的二阶初值常微分方程。<br>总结一下，二阶常微分方程的边值问题<br>$$<br>y’’ = f(x,y,y’), \quad x\in[a,b] \\<br>y(a) = \alpha, \quad y(b) = \beta<br>$$<br>的求解步骤如下：</p>
<ol>
<li>假定一个 $\gamma_1$ 值，求解初值问题<br>$$<br>y’’ = f(x,y,y’), \quad x\in[a,b] \\<br>y(a) = \alpha, \quad y’(a) = \gamma_1<br>$$<br>然后计算 $F(\gamma_1)=y(b;\gamma_1)$</li>
<li>求解关于 $W$ 的初值问题<br>$$<br>W’’=\frac{\partial f(x,y,y’)}{\partial y} W + \frac{\partial f(x,y,y’)}{\partial y’}W’ \\<br>W(a) = 0, W’(a) = 1<br>$$<br>然后计算 $F’(\gamma_1) = W(b;\gamma_1)$。</li>
<li>更新$\gamma$ 的值，<br>$$<br>\gamma_2 = \gamma_1-\frac{F(\gamma_1)}{F’(\gamma_1)}<br>$$<br>继续迭代，直到最后得到的 $y(b;\gamma_k)$ 与 $\beta$ 足够接近为止。</li>
</ol>
<h5 id="2-_OpenFOAM_中的实现">2. OpenFOAM 中的实现</h5><p>为了在OpenFOAM中求解一个任意阶常微分方程的初值问题，需要做如下准备。<br>考虑一个通用形式的常微分方程<br>$$<br>y^{n}=f(x,y,y’,\cdots,y^{n-1})<br>$$<br>定义<br>$$<br>\begin{align*}<br>y_1 &amp;=y\\<br>y_2 &amp;=y’\\<br>y_j &amp;=y^{\,j-1}, \quad j=1,2,\cdots,n<br>\end{align*}<br>$$<br>如果是求解刚性问题的 ODE 求解器，还需要定义 jacobian 矩阵。<br>令<br>$$<br>\begin{align*}<br>f_1 &amp;=y’=y_2\\<br>f_j &amp;=y’_{j}=y^{\,j+1}, \quad j=1,2,\cdots,n<br>\end{align*}<br>$$<br>则 jacobian 矩阵<br>$$<br>\begin{equation*}<br>J =<br>\begin{bmatrix}<br>\frac{\partial f_1}{\partial y_1} &amp; \frac{\partial f_1}{\partial y_2} &amp;<br>\cdots &amp;\frac{\partial f_1}{\partial y_n}\\<br>\frac{\partial f_2}{\partial y_1} &amp; \frac{\partial f_2}{\partial y_2} &amp;<br>\cdots &amp;\frac{\partial f_2}{\partial y_n}\\<br>\vdots &amp; \vdots &amp; \ddots &amp; \vdots \\<br>\frac{\partial f_n}{\partial y_1} &amp; \cdots &amp; \cdots<br>&amp;\frac{\partial f_n}{\partial y_n}<br>\end{bmatrix}<br>\end{equation*}<br>$$<br>此外，还需要给出 $f_1, f_2,\cdots,f_n$ 对自变量 $x$ 的偏导数，$\frac{\partial f_1}{\partial x}, \frac{\partial f_2}{\partial x},\cdots,\frac{\partial f_n}{\partial x}$。</p>
<p>下面举一个例子来具体说明。以常微分方程<br>$$<br>y’’=2x+2, x\in[0,1] \\<br>y(0) = 0, y’(0)= 0<br>$$<br>为例，需要定义的量为<br>$$<br>\begin{align*}<br>f_1 &amp;=y’=y_2 \\<br>f_2 &amp;= y’_{2} = 2x+2<br>\end{align*}<br>$$</p>
<p>$$<br>\begin{equation*}<br>J=<br>\begin{bmatrix}<br>0 &amp; 1 \\<br>0 &amp; 0<br>\end{bmatrix}<br>\end{equation*}<br>$$</p>
<p>$$<br>\frac{\partial f_1}{\partial x} = 0, \frac{\partial f_2}{\partial x} = 2<br>$$</p>
<p>求解这个常微分方程的代码如下：<br><figure class="highlight gherkin"><table><tr><td class="gutter"><pre><span class="line">1</span><br><span class="line">2</span><br><span class="line">3</span><br><span class="line">4</span><br><span class="line">5</span><br><span class="line">6</span><br><span class="line">7</span><br><span class="line">8</span><br><span class="line">9</span><br><span class="line">10</span><br><span class="line">11</span><br><span class="line">12</span><br><span class="line">13</span><br><span class="line">14</span><br><span class="line">15</span><br><span class="line">16</span><br><span class="line">17</span><br><span class="line">18</span><br><span class="line">19</span><br><span class="line">20</span><br><span class="line">21</span><br><span class="line">22</span><br><span class="line">23</span><br><span class="line">24</span><br><span class="line">25</span><br><span class="line">26</span><br><span class="line">27</span><br><span class="line">28</span><br><span class="line">29</span><br><span class="line">30</span><br><span class="line">31</span><br><span class="line">32</span><br><span class="line">33</span><br><span class="line">34</span><br><span class="line">35</span><br><span class="line">36</span><br><span class="line">37</span><br><span class="line">38</span><br><span class="line">39</span><br><span class="line">40</span><br><span class="line">41</span><br><span class="line">42</span><br><span class="line">43</span><br><span class="line">44</span><br><span class="line">45</span><br><span class="line">46</span><br><span class="line">47</span><br><span class="line">48</span><br><span class="line">49</span><br><span class="line">50</span><br><span class="line">51</span><br><span class="line">52</span><br><span class="line">53</span><br><span class="line">54</span><br><span class="line">55</span><br><span class="line">56</span><br><span class="line">57</span><br><span class="line">58</span><br><span class="line">59</span><br><span class="line">60</span><br><span class="line">61</span><br><span class="line">62</span><br><span class="line">63</span><br><span class="line">64</span><br><span class="line">65</span><br><span class="line">66</span><br><span class="line">67</span><br><span class="line">68</span><br><span class="line">69</span><br><span class="line">70</span><br><span class="line">71</span><br><span class="line">72</span><br><span class="line">73</span><br><span class="line">74</span><br><span class="line">75</span><br><span class="line">76</span><br><span class="line">77</span><br><span class="line">78</span><br><span class="line">79</span><br><span class="line">80</span><br><span class="line">81</span><br><span class="line">82</span><br><span class="line">83</span><br><span class="line">84</span><br><span class="line">85</span><br><span class="line">86</span><br><span class="line">87</span><br><span class="line">88</span><br><span class="line">89</span><br><span class="line">90</span><br><span class="line">91</span><br><span class="line">92</span><br><span class="line">93</span><br><span class="line">94</span><br><span class="line">95</span><br><span class="line">96</span><br><span class="line">97</span><br><span class="line">98</span><br><span class="line">99</span><br><span class="line">100</span><br><span class="line">101</span><br><span class="line">102</span><br><span class="line">103</span><br><span class="line">104</span><br><span class="line">105</span><br><span class="line">106</span><br><span class="line">107</span><br><span class="line">108</span><br><span class="line">109</span><br><span class="line">110</span><br><span class="line">111</span><br><span class="line">112</span><br><span class="line">113</span><br><span class="line">114</span><br><span class="line">115</span><br><span class="line">116</span><br><span class="line">117</span><br><span class="line">118</span><br><span class="line">119</span><br></pre></td><td class="code"><pre><span class="line">/<span class="keyword">*</span><span class="keyword">*</span><span class="keyword">*</span><span class="keyword">*</span><span class="keyword">*</span><span class="keyword">*</span><span class="keyword">*</span><span class="keyword">*</span><span class="keyword">*</span><span class="keyword">*</span><span class="keyword">*</span><span class="keyword">*</span><span class="keyword">*</span><span class="keyword">*</span><span class="keyword">*</span><span class="keyword">*</span><span class="keyword">*</span><span class="keyword">*</span><span class="keyword">*</span><span class="keyword">*</span><span class="keyword">*</span><span class="keyword">*</span><span class="keyword">*</span><span class="keyword">*</span><span class="keyword">*</span><span class="keyword">*</span><span class="keyword">*</span><span class="keyword">*</span><span class="keyword">*</span><span class="keyword">*</span><span class="keyword">*</span><span class="keyword">*</span><span class="keyword">*</span><span class="keyword">*</span><span class="keyword">*</span><span class="keyword">*</span><span class="keyword">*</span><span class="keyword">*</span><span class="keyword">*</span><span class="keyword">*</span><span class="keyword">*</span><span class="keyword">*</span><span class="keyword">*</span><span class="keyword">*</span><span class="keyword">*</span><span class="keyword">*</span><span class="keyword">*</span><span class="keyword">*</span><span class="keyword">*</span><span class="keyword">*</span><span class="keyword">*</span><span class="keyword">*</span><span class="keyword">*</span><span class="keyword">*</span><span class="keyword">*</span><span class="keyword">*</span><span class="keyword">*</span></span><br><span class="line">Description</span><br><span class="line"></span><br><span class="line">d2y/dx2 = ax + b, with a=2, b=2.</span><br><span class="line">initial value: y(0) = 0; y'(0) = 0</span><br><span class="line"></span><br><span class="line">analytical solution: y = 1/3<span class="keyword">*</span>x^3 + x^2;</span><br><span class="line"><span class="keyword">*</span><span class="keyword">*</span><span class="keyword">*</span><span class="keyword">*</span><span class="keyword">*</span><span class="keyword">*</span><span class="keyword">*</span><span class="keyword">*</span><span class="keyword">*</span><span class="keyword">*</span><span class="keyword">*</span><span class="keyword">*</span><span class="keyword">*</span><span class="keyword">*</span><span class="keyword">*</span><span class="keyword">*</span><span class="keyword">*</span><span class="keyword">*</span><span class="keyword">*</span><span class="keyword">*</span><span class="keyword">*</span><span class="keyword">*</span><span class="keyword">*</span><span class="keyword">*</span><span class="keyword">*</span><span class="keyword">*</span><span class="keyword">*</span><span class="keyword">*</span><span class="keyword">*</span><span class="keyword">*</span><span class="keyword">*</span><span class="keyword">*</span><span class="keyword">*</span><span class="keyword">*</span><span class="keyword">*</span><span class="keyword">*</span><span class="keyword">*</span><span class="keyword">*</span><span class="keyword">*</span><span class="keyword">*</span><span class="keyword">*</span><span class="keyword">*</span><span class="keyword">*</span><span class="keyword">*</span><span class="keyword">*</span><span class="keyword">*</span><span class="keyword">*</span><span class="keyword">*</span><span class="keyword">*</span><span class="keyword">*</span><span class="keyword">*</span><span class="keyword">*</span><span class="keyword">*</span><span class="keyword">*</span><span class="keyword">*</span><span class="keyword">*</span><span class="keyword">*</span><span class="keyword">*</span>/</span><br><span class="line"></span><br><span class="line"><span class="comment">#include "argList.H"</span></span><br><span class="line"><span class="comment">#include "IOmanip.H"</span></span><br><span class="line"><span class="comment">#include "ODESystem.H"</span></span><br><span class="line"><span class="comment">#include "ODESolver.H"</span></span><br><span class="line"></span><br><span class="line">using namespace Foam;</span><br><span class="line"></span><br><span class="line">·// <span class="keyword">*</span> <span class="keyword">*</span> <span class="keyword">*</span> <span class="keyword">*</span> <span class="keyword">*</span> <span class="keyword">*</span> <span class="keyword">*</span> <span class="keyword">*</span> <span class="keyword">*</span> <span class="keyword">*</span> <span class="keyword">*</span> <span class="keyword">*</span> <span class="keyword">*</span> <span class="keyword">*</span> <span class="keyword">*</span> <span class="keyword">*</span> <span class="keyword">*</span> <span class="keyword">*</span> <span class="keyword">*</span> <span class="keyword">*</span> <span class="keyword">*</span> <span class="keyword">*</span> <span class="keyword">*</span> <span class="keyword">*</span> <span class="keyword">*</span> <span class="keyword">*</span> <span class="keyword">*</span> <span class="keyword">*</span> <span class="keyword">*</span> <span class="keyword">*</span> <span class="keyword">*</span> <span class="keyword">*</span> <span class="keyword">*</span> <span class="keyword">*</span> <span class="keyword">*</span> <span class="keyword">*</span> <span class="keyword">*</span> //</span><br><span class="line"></span><br><span class="line">class myODE2</span><br><span class="line">:</span><br><span class="line">    public ODESystem</span><br><span class="line">&#123;</span><br><span class="line"></span><br><span class="line">    const scalar a_; //parameter</span><br><span class="line">    const scalar b_; //parameter</span><br><span class="line"></span><br><span class="line"></span><br><span class="line">public:</span><br><span class="line"></span><br><span class="line">    myODE2(const scalar&amp; a, const scalar&amp; b)</span><br><span class="line">	:ODESystem(),</span><br><span class="line">	a_(a),</span><br><span class="line">	b_(b)</span><br><span class="line">    &#123;&#125;</span><br><span class="line"></span><br><span class="line">    label nEqns() const // number of equations, equals to the order of ODE</span><br><span class="line">    &#123;</span><br><span class="line">        return 2;</span><br><span class="line">    &#125;</span><br><span class="line"></span><br><span class="line">    void derivatives</span><br><span class="line">    (</span><br><span class="line">        const scalar x,</span><br><span class="line">        const scalarField&amp; y,</span><br><span class="line">        scalarField&amp; dydx</span><br><span class="line">    ) const</span><br><span class="line">    &#123;</span><br><span class="line">        dydx[0] = y[1];         //f1</span><br><span class="line">        dydx[1] = a_<span class="keyword">*</span>x + b_;    //f2</span><br><span class="line">    &#125;</span><br><span class="line"></span><br><span class="line">    void jacobian // optional</span><br><span class="line">    (</span><br><span class="line">        const scalar x,</span><br><span class="line">        const scalarField&amp; y,</span><br><span class="line">        scalarField&amp; dfdx,</span><br><span class="line">        scalarSquareMatrix&amp; dfdy</span><br><span class="line">    ) const</span><br><span class="line">    &#123;</span><br><span class="line">        dfdx[0] = 0.0;        //df1/dx</span><br><span class="line">        dfdx[1] = a_;         //df2/dx</span><br><span class="line"></span><br><span class="line">        dfdy[0][0] = 0.0;     //df1/dy1</span><br><span class="line">        dfdy[0][1] = 1.0;     //df1/dy2</span><br><span class="line"></span><br><span class="line">        dfdy[1][0] = 0.0;     //df2/dy1</span><br><span class="line">        dfdy[1][1] = 0.0;     //df2/dy2</span><br><span class="line"></span><br><span class="line">    &#125;</span><br><span class="line">&#125;;</span><br><span class="line"></span><br><span class="line"></span><br><span class="line">// <span class="keyword">*</span> <span class="keyword">*</span> <span class="keyword">*</span> <span class="keyword">*</span> <span class="keyword">*</span> <span class="keyword">*</span> <span class="keyword">*</span> <span class="keyword">*</span> <span class="keyword">*</span> <span class="keyword">*</span> <span class="keyword">*</span> <span class="keyword">*</span> <span class="keyword">*</span> <span class="keyword">*</span> <span class="keyword">*</span> <span class="keyword">*</span> <span class="keyword">*</span> <span class="keyword">*</span> <span class="keyword">*</span> <span class="keyword">*</span> <span class="keyword">*</span> <span class="keyword">*</span> <span class="keyword">*</span> <span class="keyword">*</span> <span class="keyword">*</span> <span class="keyword">*</span> <span class="keyword">*</span> <span class="keyword">*</span> <span class="keyword">*</span> <span class="keyword">*</span> <span class="keyword">*</span> <span class="keyword">*</span> <span class="keyword">*</span> <span class="keyword">*</span> <span class="keyword">*</span> <span class="keyword">*</span> <span class="keyword">*</span> //</span><br><span class="line">// Main program:</span><br><span class="line"></span><br><span class="line">int main(int argc, char <span class="keyword">*</span>argv[])</span><br><span class="line">&#123;</span><br><span class="line">    argList::validArgs.append(<span class="string">"ODESolver"</span>);</span><br><span class="line">    argList args(argc, argv);</span><br><span class="line"></span><br><span class="line">    const scalar a = 2.0; </span><br><span class="line">    const scalar b = 2.0; </span><br><span class="line"></span><br><span class="line">    const label n = 100;         //number of steps</span><br><span class="line">    const scalar endTime = 1.0;  //upper bound of the interval</span><br><span class="line"></span><br><span class="line">    // Create the ODE system</span><br><span class="line">    myODE2 ode(a, b);</span><br><span class="line"></span><br><span class="line">    dictionary dict;</span><br><span class="line">    dict.add(<span class="string">"solver"</span>, args[1]);</span><br><span class="line"></span><br><span class="line">    // Create the selected ODE system solver</span><br><span class="line">    autoPtr<span class="variable">&lt;ODESolver&gt;</span> odeSolver = ODESolver::New(ode, dict);</span><br><span class="line"></span><br><span class="line">    // Initialise the ODE system fields</span><br><span class="line">    scalar xStart = 0.0;        // lower bound of the interval</span><br><span class="line">    scalar dx = endTime/n;      //step value</span><br><span class="line"></span><br><span class="line">    scalarField yStart(ode.nEqns());</span><br><span class="line">    yStart[0] = 0.0; // initial value of y</span><br><span class="line">    yStart[1] = 0.0; // initial value of y'</span><br><span class="line"></span><br><span class="line">    scalar dxEst = 0.1;</span><br><span class="line">    scalar xEnd  = 0.0;</span><br><span class="line"></span><br><span class="line">    scalarField dyStart(ode.nEqns()); // dyStart[0]=f1, dyStart[1]=f2 ...</span><br><span class="line"></span><br><span class="line">    for(label i =0; i<span class="variable">&lt;n; i++)</span><br><span class="line">    &#123;</span><br><span class="line">	    xEnd = xStart + dx;</span><br><span class="line">	    ode.derivatives(xStart, yStart, dyStart);</span><br><span class="line">	    odeSolver-&gt;</span>solve(xStart, xEnd, yStart, dxEst);</span><br><span class="line">	    xStart = xEnd;</span><br><span class="line">	    Info <span class="variable">&lt;&lt; xStart &lt;&lt; "    " &lt;&lt; yStart[0] &lt;&lt; endl; // output (x,y) for each dx.</span><br><span class="line">    &#125;</span><br><span class="line"></span><br><span class="line">    return 0;</span><br><span class="line">&#125;</span></span><br></pre></td></tr></table></figure></p>
<p>编译之后，假设你的可执行程序名为 <code>TestODE</code>，则运行<br><figure class="highlight applescript"><table><tr><td class="gutter"><pre><span class="line">1</span><br></pre></td><td class="code"><pre><span class="line">TestODE RKCK45 &gt; <span class="command">log</span></span><br></pre></td></tr></table></figure></p>
<p>就得到了数值解。<br>将数值解与解析解画图如下<br><img src="/image/RKCK45.png" alt=""><br>可见在这里简单例子中，数值解与解析解吻合非常好。<br>同时注意，这个例子，用显式欧拉方法无法得到收敛的解。</p>
<p><strong>参考资料</strong>：</p>
<ol>
<li><a href="http://hassankassem.me/posts/ode/" target="_blank" rel="external">http://hassankassem.me/posts/ode/</a></li>
<li>余德浩 , 汤华中， 微分方程数值解法，科学出版社，2003</li>
</ol>

      
    </div>
    <footer class="article-footer">
      <a data-url="http://xiaopingqiu.github.io/2016/08/21/ODE/" data-id="cj2lay37b0062kkmb2aup5a42" class="article-share-link">Share</a>
      
      
  <ul class="article-tag-list"><li class="article-tag-list-item"><a class="article-tag-list-link" href="/tags/ODE/">ODE</a></li></ul>

    </footer>
  </div>
  
    
<nav id="article-nav">
  
    <a href="/2016/08/27/ParaviewStreamLineOnSlice/" id="article-nav-newer" class="article-nav-link-wrap">
      <strong class="article-nav-caption">Newer</strong>
      <div class="article-nav-title">
        
          在 Paraview 中画截面上的流线
        
      </div>
    </a>
  
  
    <a href="/2016/08/21/fixedFluxPressure/" id="article-nav-older" class="article-nav-link-wrap">
      <strong class="article-nav-caption">Older</strong>
      <div class="article-nav-title">fixedFluxPressure 边界条件</div>
    </a>
  
</nav>

  
</article>


<section id="comments">
  <!-- 多说评论框 start -->
  <div class="ds-thread" data-thread-key="post-ODE" data-title="OpenFOAM 中求解 ODE 一例" data-url="http://xiaopingqiu.github.io/2016/08/21/ODE/"></div>
  <!-- 多说评论框 end -->
  <!-- 多说公共JS代码 start (一个网页只需插入一次) -->
  <script type="text/javascript">
  var duoshuoQuery = {short_name:"xiaopingqiu"};
	(function() {
		var ds = document.createElement('script');
		ds.type = 'text/javascript';ds.async = true;
		ds.src = (document.location.protocol == 'https:' ? 'https:' : 'http:') + '//static.duoshuo.com/embed.js';
		ds.charset = 'UTF-8';
		(document.getElementsByTagName('head')[0] 
		 || document.getElementsByTagName('body')[0]).appendChild(ds);
	})();
	</script>
  <!-- 多说公共JS代码 end -->
</section>

</section>
        
          <aside id="sidebar">
  
    
  <div class="widget-wrap">
    <h3 class="widget-title">Categories</h3>
    <div class="widget">
      <ul class="category-list"><li class="category-list-item"><a class="category-list-link" href="/categories/C/">C++</a><span class="category-list-count">2</span></li><li class="category-list-item"><a class="category-list-link" href="/categories/DEM/">DEM</a><span class="category-list-count">1</span></li><li class="category-list-item"><a class="category-list-link" href="/categories/OpenFOAM/">OpenFOAM</a><span class="category-list-count">44</span></li><li class="category-list-item"><a class="category-list-link" href="/categories/Paraview/">Paraview</a><span class="category-list-count">5</span></li><li class="category-list-item"><a class="category-list-link" href="/categories/swak4Foam/">swak4Foam</a><span class="category-list-count">1</span></li><li class="category-list-item"><a class="category-list-link" href="/categories/test/">test</a><span class="category-list-count">2</span></li><li class="category-list-item"><a class="category-list-link" href="/categories/vim/">vim</a><span class="category-list-count">1</span></li></ul>
    </div>
  </div>

  
    
  <div class="widget-wrap">
    <h3 class="widget-title">Tags</h3>
    <div class="widget">
      <ul class="tag-list"><li class="tag-list-item"><a class="tag-list-link" href="/tags/Boundary-conditions/">Boundary conditions</a><span class="tag-list-count">6</span></li><li class="tag-list-item"><a class="tag-list-link" href="/tags/C/">C++</a><span class="tag-list-count">2</span></li><li class="tag-list-item"><a class="tag-list-link" href="/tags/CentOS/">CentOS</a><span class="tag-list-count">1</span></li><li class="tag-list-item"><a class="tag-list-link" href="/tags/Code-Explained/">Code Explained</a><span class="tag-list-count">29</span></li><li class="tag-list-item"><a class="tag-list-link" href="/tags/LES/">LES</a><span class="tag-list-count">1</span></li><li class="tag-list-item"><a class="tag-list-link" href="/tags/LIGGGHTS/">LIGGGHTS</a><span class="tag-list-count">1</span></li><li class="tag-list-item"><a class="tag-list-link" href="/tags/ODE/">ODE</a><span class="tag-list-count">1</span></li><li class="tag-list-item"><a class="tag-list-link" href="/tags/OpenFOAM/">OpenFOAM</a><span class="tag-list-count">20</span></li><li class="tag-list-item"><a class="tag-list-link" href="/tags/Postprocessing/">Postprocessing</a><span class="tag-list-count">9</span></li><li class="tag-list-item"><a class="tag-list-link" href="/tags/Preprocessing/">Preprocessing</a><span class="tag-list-count">2</span></li><li class="tag-list-item"><a class="tag-list-link" href="/tags/RTS/">RTS</a><span class="tag-list-count">3</span></li><li class="tag-list-item"><a class="tag-list-link" href="/tags/TIL/">TIL</a><span class="tag-list-count">1</span></li><li class="tag-list-item"><a class="tag-list-link" href="/tags/Windows/">Windows</a><span class="tag-list-count">1</span></li><li class="tag-list-item"><a class="tag-list-link" href="/tags/fvOptions/">fvOptions</a><span class="tag-list-count">2</span></li><li class="tag-list-item"><a class="tag-list-link" href="/tags/groovyBC/">groovyBC</a><span class="tag-list-count">1</span></li><li class="tag-list-item"><a class="tag-list-link" href="/tags/paraview/">paraview</a><span class="tag-list-count">1</span></li><li class="tag-list-item"><a class="tag-list-link" href="/tags/test/">test</a><span class="tag-list-count">2</span></li><li class="tag-list-item"><a class="tag-list-link" href="/tags/thermophysicalModels/">thermophysicalModels</a><span class="tag-list-count">5</span></li><li class="tag-list-item"><a class="tag-list-link" href="/tags/turbulence-model/">turbulence model</a><span class="tag-list-count">7</span></li><li class="tag-list-item"><a class="tag-list-link" href="/tags/vim/">vim</a><span class="tag-list-count">1</span></li><li class="tag-list-item"><a class="tag-list-link" href="/tags/wall-functions/">wall functions</a><span class="tag-list-count">4</span></li></ul>
    </div>
  </div>

  
    
  <div class="widget-wrap">
    <h3 class="widget-title">Tag Cloud</h3>
    <div class="widget tagcloud">
      <a href="/tags/Boundary-conditions/" style="font-size: 15.56px;">Boundary conditions</a><a href="/tags/C/" style="font-size: 11.11px;">C++</a><a href="/tags/CentOS/" style="font-size: 10px;">CentOS</a><a href="/tags/Code-Explained/" style="font-size: 20px;">Code Explained</a><a href="/tags/LES/" style="font-size: 10px;">LES</a><a href="/tags/LIGGGHTS/" style="font-size: 10px;">LIGGGHTS</a><a href="/tags/ODE/" style="font-size: 10px;">ODE</a><a href="/tags/OpenFOAM/" style="font-size: 18.89px;">OpenFOAM</a><a href="/tags/Postprocessing/" style="font-size: 17.78px;">Postprocessing</a><a href="/tags/Preprocessing/" style="font-size: 11.11px;">Preprocessing</a><a href="/tags/RTS/" style="font-size: 12.22px;">RTS</a><a href="/tags/TIL/" style="font-size: 10px;">TIL</a><a href="/tags/Windows/" style="font-size: 10px;">Windows</a><a href="/tags/fvOptions/" style="font-size: 11.11px;">fvOptions</a><a href="/tags/groovyBC/" style="font-size: 10px;">groovyBC</a><a href="/tags/paraview/" style="font-size: 10px;">paraview</a><a href="/tags/test/" style="font-size: 11.11px;">test</a><a href="/tags/thermophysicalModels/" style="font-size: 14.44px;">thermophysicalModels</a><a href="/tags/turbulence-model/" style="font-size: 16.67px;">turbulence model</a><a href="/tags/vim/" style="font-size: 10px;">vim</a><a href="/tags/wall-functions/" style="font-size: 13.33px;">wall functions</a>
    </div>
  </div>

  
    
  <div class="widget-wrap">
    <h3 class="widget-title">Archives</h3>
    <div class="widget">
      <ul class="archive-list"><li class="archive-list-item"><a class="archive-list-link" href="/archives/2017/05/">五月 2017</a><span class="archive-list-count">1</span></li><li class="archive-list-item"><a class="archive-list-link" href="/archives/2016/08/">八月 2016</a><span class="archive-list-count">8</span></li><li class="archive-list-item"><a class="archive-list-link" href="/archives/2016/06/">六月 2016</a><span class="archive-list-count">5</span></li><li class="archive-list-item"><a class="archive-list-link" href="/archives/2016/05/">五月 2016</a><span class="archive-list-count">3</span></li><li class="archive-list-item"><a class="archive-list-link" href="/archives/2016/04/">四月 2016</a><span class="archive-list-count">12</span></li><li class="archive-list-item"><a class="archive-list-link" href="/archives/2016/03/">三月 2016</a><span class="archive-list-count">6</span></li><li class="archive-list-item"><a class="archive-list-link" href="/archives/2015/12/">十二月 2015</a><span class="archive-list-count">1</span></li><li class="archive-list-item"><a class="archive-list-link" href="/archives/2015/11/">十一月 2015</a><span class="archive-list-count">2</span></li><li class="archive-list-item"><a class="archive-list-link" href="/archives/2015/10/">十月 2015</a><span class="archive-list-count">1</span></li><li class="archive-list-item"><a class="archive-list-link" href="/archives/2015/09/">九月 2015</a><span class="archive-list-count">6</span></li><li class="archive-list-item"><a class="archive-list-link" href="/archives/2015/08/">八月 2015</a><span class="archive-list-count">1</span></li><li class="archive-list-item"><a class="archive-list-link" href="/archives/2015/06/">六月 2015</a><span class="archive-list-count">2</span></li><li class="archive-list-item"><a class="archive-list-link" href="/archives/2015/05/">五月 2015</a><span class="archive-list-count">6</span></li><li class="archive-list-item"><a class="archive-list-link" href="/archives/2015/04/">四月 2015</a><span class="archive-list-count">2</span></li></ul>
    </div>
  </div>

  
    
  <div class="widget-wrap">
    <h3 class="widget-title">Recents</h3>
    <div class="widget">
      <ul>
        
          <li>
            <a href="/2017/05/12/duoshuoAnnouncement/">多说评论系统将停止提供服务</a>
          </li>
        
          <li>
            <a href="/2016/08/27/ParaviewScritps/">Paraview 脚本一例</a>
          </li>
        
          <li>
            <a href="/2016/08/27/ParaviewCamera/">Paraview 中有关 Camera 的操作两例</a>
          </li>
        
          <li>
            <a href="/2016/08/27/ParaviewCustomFilter/">Paraview 中创建 Custom Filter</a>
          </li>
        
          <li>
            <a href="/2016/08/27/ParaviewStreamLineOnSlice/">在 Paraview 中画截面上的流线</a>
          </li>
        
      </ul>
    </div>
  </div>

  
</aside>
        
      </div>
      <footer id="footer">
  
  <div class="outer">
    <div id="footer-info" class="inner">
      &copy; 2017 Giskard Q.<br>
      Powered by <a href="http://hexo.io/" target="_blank">Hexo</a>
    </div>
  </div>
</footer>

<script src="//dn-lbstatics.qbox.me/lbservice/busuanzi/2.0/busuanzi.mini.js"/></script>

    </div>
    <nav id="mobile-nav">
  
    <a href="/" class="mobile-nav-link">Home</a>
  
    <a href="/archives" class="mobile-nav-link">Archives</a>
  
    <a href="/atom.xml" class="mobile-nav-link">Rss</a>
  
    <a href="/about" class="mobile-nav-link">About</a>
  
</nav>
    

<script src="//ajax.googleapis.com/ajax/libs/jquery/2.0.3/jquery.min.js"></script>


  <link rel="stylesheet" href="/fancybox/jquery.fancybox.css" type="text/css">
  <script src="/fancybox/jquery.fancybox.pack.js" type="text/javascript"></script>


<script src="/js/script.js" type="text/javascript"></script>

<script>
  (function(i,s,o,g,r,a,m){i['GoogleAnalyticsObject']=r;i[r]=i[r]||function(){
  (i[r].q=i[r].q||[]).push(arguments)},i[r].l=1*new Date();a=s.createElement(o),
  m=s.getElementsByTagName(o)[0];a.async=1;a.src=g;m.parentNode.insertBefore(a,m)
  })(window,document,'script','//www.google-analytics.com/analytics.js','ga');

  ga('create', 'UA-62501686-1', 'auto');
  ga('send', 'pageview');

</script>

  </div>
<!-- mathjax config similar to math.stackexchange -->

<script type="text/x-mathjax-config">
  MathJax.Hub.Config({
    tex2jax: {
      inlineMath: [ ['$','$'], ["\\(","\\)"] ],
      processEscapes: true
    }
  });
</script>

<script type="text/x-mathjax-config">
    MathJax.Hub.Config({
      tex2jax: {
        skipTags: ['script', 'noscript', 'style', 'textarea', 'pre', 'code']
      }
    });
</script>

<script type="text/x-mathjax-config">
    MathJax.Hub.Queue(function() {
        var all = MathJax.Hub.getAllJax(), i;
        for(i=0; i < all.length; i += 1) {
            all[i].SourceElement().parentNode.className += ' has-jax';
        }
    });
</script>

<script type="text/javascript" src="http://cdn.mathjax.org/mathjax/latest/MathJax.js?config=TeX-AMS-MML_HTMLorMML">
</script>
</body>
</html>